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In the framework of a recently reported linear-scaling method for density-functional- 
pseudopotential calculations, we investigate the use of localized basis functions for such work. We 
propose a basis set in which each local orbital is represented in terms of an array of "blip functions" 
on the points of a grid. We analyze the relation between blip-function basis sets and the plane-wave 
basis used in standard pseudopotential methods, derive criteria for the approximate equivalence 
of the two, and describe practical tests of these criteria. Techniques are presented for using blip- 
function basis sets in linear-scaling calculations, and numerical tests of these techniques are reported 
for Si crystal using both local and non-local pseudopotentials. We find rapid convergence of the total 
energy to the values given by standard plane-wave calculations as the radius of the linear-scaling 
localized orbitals is increased. 



1 



I. INTRODUCTION 

First-principles calculations based on density func- 
tional theorjocl (DFT) and the pseudopotential method 
are widely used for studying tho-energetics, structure and 
dynamics of solids and liquidstTtl. In the standard ap- 
proach, the occupied Kohn-Sham orbitals are expanded 
in terms of plane waves, and the ground state is found 
by minimizing the total energy with respect to the plane- 
wave coefficients^]. Calculations on systems of over a 
hundred atoms with this approach are now quite com- 
mon. However, it has proved difficult to go to very much 
larger systems, because the computational effort in this 
approach depends on the number of atoms N at least 
as A 2 , and asymptotically as N 3 . Because of this limi- 
tation, there has been a vigorous effort jn— the past few 
years to develop linear-scaling methods^! 2 -! - methods 
in which the effort depends only linearly on the number 
of atoms. We have recently described a general theoreti- 
cal frameworkibjc-developing linear-scaling self-consistent 
DFT schemesEna. We presented one practical way of 
implementing such a scheme, and investigated its perfor- 
mance for crystalline silicon. Closely related ideas have 
been reported by other authorsExcil - an overview of work 
on linear-scaling methods was given in the Introduction 
of our previous papeiH. 

The practical feasibility of linear-scaling DFT tech- 
niques is thus well established. However, there are still 
technical problems to be solved before the techniques can 
be routinely applied. Our aim here is to study the prob- 
lem of representing the localized orbitals that appear in 
linear-scaling methods (support functions in our termi- 
nology) - in other words, the problem of basis functions. 
To put this in context, we recall briefly the main ideas of 
our linear-scaling DFT method. 

Standard DFT can be expressed in terms of the Kohn- 
Sham density matrix p(r, r') . The total-energy functional 
can be written in terms of p, and the ground state is 
obtained by minimization with respect to p subject to 
two constraints: p is idempotent (it is a projector, so 
that its eigenvalues are or 1), and its trace is equal to 
half the number of electrons. Linear-scaling behavior is 
obtained by imposing a limitation on the spatial range of 
P- 

p(r,r') = 0, \v-r'\>R c . (1) 

By the variational principle, we then get an upper bound 
E(R C ) to the true ground-state energy E . Since the true 
ground-state density matrix decays to zero as r — r' | — > 
oo, we expect that E(R C — > oo) = Eq. To make the 
scheme practicable, we introduced the further condition 
that p be separable: 

p(ry) = J2ct> a (r)K a p<f> (v') , (2) 

a/3 

where the number of support functions Q (r) is finite. 
The limitation on the spatial range of p is imposed by 



requiring that the <fi a (r) are non-zero only in localized 
regions ( "support regions" ) and that the spatial range of 
K a p is limited. In our method, the support regions are 
centered on the at oms , and move with them. 

We have shownEiEj that the condition on the eigen- 
values of p canJae satisfied by the method of Li, Nunes 
and VanderbiltO (LNV): instead of directly varying p, 
we express it as: 

p = 3cr * a — 2a * a * a , (3) 

where the asterisk indicates the continuum analog of ma- 
trix multiplication. As shown by LNV, this representa- 
tion of p not only ensures that its eigenvalues lie in the 
range [0, 1], but it drives them towards the values and 
1. In our scheme, the auxiliary matrix cr(r, r') has the 
same type of separability as p: 

< 7(r,r')=^^(r)i^^(r'). (4) 

a/3 

This means that K is given by the matrix equation: 

K = 3LSL - 2LSLSL , (5) 
where Sap is the overlap matrix of support functions: 

S a p = fdr<f> a <pp . (6) 

We can therefore summarize the overall scheme as fol- 
lows. The total energy is expressed in terms of p, which 
depends on the separable quantity a. The ground-state 
energy is obtained by minimization with respect to the 
support functions Q (r) and the matrix elements L a p, 
with the 4> a confined to localized regions centered on the 
atoms, and the L a p subject to a spatial cut-off. 

The ^a(r) must be allowed to vary freely in the min- 
imization process, just like the Kohn-Sham orbitals in 
conventional DFT, and we must consider how to repre- 
sent them. As always, there is a choice: we pGap-rep- 
resent them either by their values on a gridcTEJ, or 
in terms .of, some set of basis functions. In our previ- 
ous workc^ES, we used a grid representation. This was 
satisfactory for discussing the feasibility of linear-scaling 
schemes, but seems to us to suffer from significant draw- 
backs. Since the support regions are centered on the 
ions in our method, this means that when the ions move, 
the boundaries of the regions will cross the grid points. 
In any simple grid-based scheme, this will cause trou- 
blesome discontinuities. In addition, the finite-difference 
representation of the kinetic-energy operator in a grid 
representation causes problems at the boundaries of the 
regions. A further point is that in a purely grid-based 
method we are almost certainly using more variables than 
are really necessary. These problems have led us to con- 
sider basis-function methods. 

We describe in this paper a practical basis-function 
scheme for linear-scaling DFT, and we study its perfor- 
mance in numerical calculations. The basis consists of 
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an array of localized functions - we call them "blip func- 
tions" . There is an array of blip functions for each sup- 
port region, and the array moves with the region. The 
use of such arrays of localized functions as a basis for 
quantum calculations is not newc3~L3. However, to our 
knowledge it has not been discussed before in the context 
of linear-scaling calculations. 

The plan of the paper is as follows. In Sec. 2, we 
emphasize the importance of considering the relation be- 
tween blip- function and plane- wave basis sets, and we use 
this relation to analyze how the calculated ground-state 
energy will depend on the width and spacing of the blip 
functions. We note some advantages of using B-splines as 
blip functions, and we then present some practical tests 
which illustrate the convergence of the ground-state en- 
ergy with respect to blip width and spacing. We then 
go on (Sec. 3) to discuss the technical problems of using 
blip-function basis sets in linear-scaling DFT. We report 
the results of practical tests, which show explicitly how 
the ground-state energy in linear-scaling DFT converges 
to the value obtained in a standard plane-wave calcu- 
lation. Section 4 gives a discussion of the results, and 
presents our conclusions. Some mathematical derivations 
are given in an appendix. 

II. BLIP FUNCTIONS AND PLANE WAVES 
A. General considerations 

Before we focus on linear-scaling problems, we need to 
set down some elementary ideas about basis functions. 
To start with, we therefore ignore the linear-scaling as- 
pects, and we discuss the general problem of solving 
Schrodinger's equation using basis functions. It is enough 
to discuss this in one dimension, and we assume a period- 
ically repeating system, so that the potential V(x) acting 
on the electrons is periodic: V(x + t) = V(x), where t 
is any translation vector. Self-consistency questions are 
irrelevant at this stage, so that V(x) is given. The gener- 
alization to three-dimensional self-consistent calculations 
will be straightforward. 

In a plane- wave basis, the wavefunctions tpi{x) are ex- 
panded as 

rl>i{x) = L- 1 ' 2 c,g exp(zGx) , (7) 

G 

where the reciprocal lattice vectors of the repeating ge- 
ometry are given by G — 2nn/L (n is an integer), and 
we include all G up to some cut-off G max . We obtain 
the ground-state energy E(G max ) in the given basis by 
minimization with respect to the qg, subject to the con- 
straints of orthonormality. For the usual variational rea- 
sons, E(G max ) is a monotonically decreasing function of 
G m ax which tends to the exact value Eq as G max — > oo. 

Now instead of plane waves we want to use an array 
of spatially localized basis functions (blip functions). Let 



fo(x) be some localized function, and denote by fi{x) the 
translated function f (x — la) , where I is an integer and a 
is a spacing, which is chosen so that L is an exact multiple 
of a: L — Ma. We use the array of blip functions fi(x) as 
a basis. Equivalently, we can use any independent linear 
combinations of the fe{x). In considering the relation 
between blip functions and plane waves, it is particularly 
convenient to work with "blip waves", xg(x), defined as: 

M-l 

X g(x) =A g J2 M x ) exp{iGRi) , (8) 

where Ri = la, and Aq is some normalization constant. 

The relation between blip waves and plane waves can 
be analyzed by considering the Fourier representation of 
Xg{x). It is straightforward to show that Xg(%) has 
Fourier components only at wavevectors G + T, where T 
is a reciprocal lattice vector of the blip grid: T = 2irm/a 
(m is an integer). In fact: 

Xg(x) = (A G /a) f(G + r) exp (i(G + I» , (9) 
r 

where f{q) is the Fourier transform of fo(x): 

/oo 
dxf (x)e^. (10) 
-oo 

At this point, it is useful to note that for some choices 
of f (x) the blip-function basis set is exactly equivalent to 
a plane-wave basis set. For this to happen, f(q) must be 
exactly zero beyond some cut-off wavevector q cu t- Then 
provided g cut > G max and provided q cut + G max < 27r/a, 
all the r / terms in Eq. (|J) will vanish and all blip 
waves for — G max < G < G max will be identical to plane- 
waves. (Of course, we must also require that f(q) ^ 
for \q\ < G max .) 

Our main aim in this Section is to determine how the 
total energy converges to the exact value as the width 
and spacing of the blip functions are varied. The spacing 
is controlled by varying a, and the width is controlled by 
scaling each blip function: fo(x) — > fo(sx), where s is 
a scaling factor. In the case of blip functions for which 
f(q) cuts off in the way just described, the convergence 
of the total energy is easy to describe. Suppose we take 
a fixed blip width, and hence a fixed wavevector cut-off 
q cut . If the blip spacing a is small enough so that q cut < 
7T /a, then it follows from what we have said that the blip 
basis set is exactly equivalent to a plane-wave basis set 
having G max = q cu t- This means that the total energy is 
equal to E(q cut ) and is completely independent of a when 
the latter falls below the threshold value ath = ^/icnt- 
This is connected with the fact that the blip basis set 
becomes over-complete when a < a t h- there are linear 
dependences between the M blip functions fe(x). 

It follows from this that the behavior of the total en- 
ergy as a function of blip spacing and blip width is as 
shown schematically in Fig. 1. As the width is reduced, 
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the cut-off g C ut increases in proportion to the scaling fac- 
tor s, so that the threshold spacing ath is proportional 
to the width. The energy value E(q cut ) obtained for 
a < ath decreases monotonically with the width, as fol- 
lows from the monotonic decrease of -E(G max ) with G max 
for a plane- wave basis set. Note that in Fig. 1 we have 
shown E at fixed width as decreasing monotonically with 
a for a > a t h- In fact, this may not always happen. De- 
crease of a does not correspond simply to addition of 
basis functions and hence to increase of variational free- 
dom: it also involves relocation of the basis functions. 
However, what is true is that E for a > ath is always 
greater than E for a < a t h, as can be proved from the 
over-completeness of the blip basis set for a < a t h- At 
large spacings, the large-width blip basis is expected to 
give the lower energy, since in this region the poorer rep- 
resentation of long waves should be the dominant source 
of error. 
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FIG. 1. Expected schematic form for the total ground-state 
energy as a function of the blip-grid spacing for two different 
blip widths, in the case where the Fourier components of the 
blip functions vanish beyond some cut-off. The horizontal 
dotted line shows the exact ground-state energy Eq. The 
vertical dashed lines mark the threshold values a of blip grid 
spacing (see text). 

Up to now, we have considered only the rather artificial 
case where the Fourier components of the blip function 
are strictly zero beyond a cut-off. This means that the 
blip function must extend to infinity in real space, and 
this is clearly no use if we wish to do all calculations in 
real space. We actually want fo(x) to be strictly zero 
beyond some real-space cut-off b a : fo(x) = if |x| > bo. 
This means that f(q) will extend to infinity in reciprocal 
space. However, we can expect that with a judicious 
choice for the form of fo{x) the Fourier components /(<?) 
will still fall off very rapidly, so that the behavior of the 
total energy is still essentially as shown in Fig. 1. If 
the choice is not judicious, we shall need a considerably 
greater effort to bring E within a specified tolerance of 
the exact value than if we were using a plane-wave basis 
set. With a plane-wave basis, a certain cut-off G max is 
needed in order to achieve a specified tolerance in E. 



With blip functions whose Fourier components cut off at 
G max , we should need a blip spacing of 7r/G max to achieve 
the same tolerance. Our requirement on the actual choice 
of blip function is that the spacing needed to achieve the 
given tolerance should be not much less than 7r/G max . 

B. B-splines as blip functions 

Given that the blip function cuts off in real space at 
some distance bo, it is helpful if the function and some 
of its derivatives go smoothly to zero at this distance. If 
fo(x) and all its derivatives up to and including the nth 
vanish at \x\ = bo, then f(q) falls off asymptotically as 
l/|g| n+2 as \q\ — > oo. One way of making a given set 
of derivatives vanish is to build fo{x) piecewise out of 
suitable polynomials. As an example, we examine here 
the choice of fo(x) as a B-spline. 

B-splines are localized polynomial basis functions that 
are equivalent to a representation of functions in terms 
of cubic splines. A single B-spline B(x) centered at the 
origin and covering the region \x\ < 2 is built out of third- 
degree polynomials in the four intervals —2 < x < — 1, 
— 1<£<0, 0<£<1 and 1 < x < 2, and is defined as: 

( 1 - |.t 2 + ||x| 3 if < |a;| < 1 
B(x) = { \{2-\x\f if 1<H<2 (11) 
[ if 2 < \x\ 

The function and its first two derivatives are continuous 
everywhere. The Fourier transform of B(x), defined as 
in Eq. @, is: 

B(q) = \(3-4cosq + cos2q) , (12) 

q 

which falls off asymptotically as g 4 , as expected. Our 
choice of blip function is thus fo(x) = B(2x/bo), so that 
f(q) = B(±b q). 

The transform B(q) falls rapidly to small values in the 
region \q\ ~ 7r. It is exactly zero at the set of wavevectors 
q n = 2im (n is a non-zero integer), and is very small 
in a rather broad region around each q n , because the 
lowest non-vanishing term in a polynomial expansion of 
3 — 4 cos q + cos 2q is of degree q 4 . This suggests that this 
choice of blip function will behave rather similarly to one 
having a Fourier cut-off q cu t = 2it /bo- In other words, 
if we keep bo fixed and reduce the blip spacing a, the 
energy should approach the value obtained in a plane- 
wave calculation having G max = g cu t when a ~ \i>o- The 
practical tests that now follow will confirm this. (We note 
that B-splines are usually employed with a blip spacing 
equal to i&o! here, however, we are allowing the spacing 
to vary). 
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C. Practical tests 

Up to now, it was convenient to work in one dimension, 
but for practical tests we clearly want to go to real three- 
dimensional systems. To do this, we simply take the blip 
function /o(r) to be the product of factors depending on 
the three Cartesian components of r: 

/ (r) = p {x) p (y) p (z) . (13) 

All the considerations outlined above for a blip basis in 
one dimension apply unchanged to the individual factors 
Po(x) etc, which are taken here to be B-splines. Corre- 
sponding to the blip grid of spacing a in one dimension, 
the blip functions /o(r) now sit on the points of a three- 
dimensional grid which we assume here to be simple cu- 
bic. The statements made above about the properties 
of the B-spline basis are expected to remain true in this 
three-dimensional form. 

We present here some tests on the performance of a 
B-spline basis for crystalline Si. At first sight, it might 
appear necessary to write a new code in order to per- 
form such tests. However, it turns out that rather minor 
modifications to an existing plane-wave code allow one 
to produce results that are identical to those that would 
be obtained with a B-spline basis. For the purpose of the 
present tests, this is sufficient. The notion behind this 
device is that blip functions can be expanded in terms 
of plane waves, so that the function-space spanned by 
a blip basis is contained within the space spanned by a 
plane-wave basis, provided the latter has a large enough 
G max . Then all we have to do to get the blip basis is 
to project from the large plane-wave space into the blip 
space. Mathematical details of how to do this projection 
in practice are given in the Appendix. The practical tests 
have been done with the CASTEP codeQ, which we have 
modified to perform the necessary projections. 

Our tests have been done on the diamond-structure 
Si crystal, using the Appelbaum-Hamanncfj local pseu- 
dopotential; this is an empirical pseudopotential, but suf- 
fices to illustrate the points of principle at issue here. 
The choice of fc-point sampling is not expected to make 
much difference to the performance of blip functions, 
and we have done the calculations with a /c-point set 
corresponding to the lowest-order 4 fc-point Monkhorst- 
Packcil sampling for an 8-atom cubic cell. If we go the 
next-order set of 32 /c-points, the total energy per atom 
changes by less than 0.1 eV. For reference purposes, we 
have first used CASTEP in its normal unmodified plane- 
wave form to examine the convergence of total energy 
with respect to plane-wave cut-off for the Appelbaum- 
Hamann potential. We find that for plane-wave cut- 
off energies E pw = ?i 2 G^ lax /2m equal to 150, 250 and 
350 eV, the total energies per atom are -115.52, -115.64 
and -115.65 eV. This means that to obtain an accuracy 
of 0.1 eV/atom, a cut-off of 150 eV (corresponding to 
Gmax = 6.31 A -1 ) is adequate. According to the discus- 
sion of Sec. 2.2, the properties of this plane-wave basis 



should be quite well reproduced by a blip-function basis 
of B-splines having half- width bo = 27r/G max = 1-0 A, 
and the total energy calculated with this basis should 
converge rapidly to the plane-wave result when the blip 
spacing falls below a ~ ^bo — 0.5 A. 
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FIG. 2. Convergence of the total energy for two different 
blip half- widths as a function of the blip grid spacing. The cal- 
culations were performed with the Appelbaum-Hamann pseu- 
dopotential. E(plw) is the plane-wave result, obtained with 
a cutoff of 250 eV. 

We have done the tests with two widths of B-splines: 
b a = 1.25 and 1.0 A, and in each case we have calculated 
E as a function of the blip spacing a. In all cases, we have 
used a plane- wave cut-off large enough to ensure that er- 
rors in representing the blip-function basis are negligible 
compared withe errors attributable to the basis itself. 
Our results for £ as a function of a (Fig. 2) fully con- 
firm our expectations. First, they have the general form 
indicated in Fig. 1. The difference is that since we have 
no sharp cut-off in reciprocal space, E does not become 
constant when a falls below a threshold, but instead con- 
tinues to decrease towards the exact value. Second, for 
b = 1.0 A, E does indeed converge rapidly to the plane- 
wave result when a falls below ~ 0.5 A. Third, the larger 
blip width gives the lower energy at larger spacings. 

III. THE BLIP-FUNCTION BASIS IN 
LINEAR-SCALING CALCULATIONS 

A. Technical matters 

In our linear-scaling scheme, the support functions 
4> a (r) must be varied within support regions, which are 
centered on the atoms. These regions, taken to be spher- 
ical with radius i? reg in our previous work, move with 
the atoms. In the present work, the (j> a are represented 
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in terms of blip functions. Each atom has a blip grid at- 
tached to it, and this grid moves rigidly with the atom. 
The blip functions sit on the point of this moving grid. 
To make the region localized, the set of blip-grid points 
is restricted to those for which the associated blip func- 
tion is wholly contained within the region radius i? r cg- If 
we denote by f a e,( r ) the £th blip function in the region 
supporting (j) a , then the representation is: 



<Mr) 



b a ef a e(r) 



(14) 



and the blip coefficients b a i have to be varied to minimize 
the total energy. 

The <p a enter the calculation through their overlap ma- 
trix elements and the matrix elements of kinetic and po- 
tential energies. The overlap matrix S a p [see Eq. (^)] can 
be expressed analytically in terms of the blip coefficients: 



S a = 2. bat bf3£> Salfigj 



(15) 



where s a t,j3i> is the overlap matrix between blip func- 
tions: 



(16) 



SalM' — / dr f a t fp£i , 



which is known analytically. Similarly, the kinetic energy 
matrix elements: 



(17) 



can be calculated analytically by writing: 

T a p = b a( bpi> t a £jjfj , (18) 



where: 



tn 



2m 



dr/ a ,V 2 



(19) 



However, matrix elements of the potential energy can- 
not be treated analytically, and their integrations must 
be approximated by summation on a grid. This 'integra- 
tion grid' is, of course, completely distinct from the blip 
grids. It does not move with the atoms, but is a single 
grid fixed in space. If the position of the mth point on 
the integration grid is called r m , then the matrix elements 
of the local potential (pseudopotential plus Hartree and 
exchange-correlation potentials) are approximated by: 

V a () = dr (p a V(pp ~ SUiat ^ <j>a(r m )V(T m )<f)p(T m ) , 

(20) 

where Su>- ln t is the volume per grid point. For a non-local 
pseudopotential, we assume the real-space version of the 



Kleinman-BylanderH representation, and the terms in 
this are also calculated as a sum over points of the inte- 
gration grid. We note that the approximate equivalence 
of the B-spline and plane-wave bases discussed above 
gives us an expectation for the required integration-grid 
spacing fe. In a plane- wave calculation, h should in prin- 
ciple be less than 7r(2G max ). But the blip-grid spacing is 
approximately a — 7rG max . We therefore expect to need 
h ks ha. 

In order to calculate V a p like this, we have to know all 
the values 4> a {r m ) on the integration grid: 



<t>a{*m) 



b a t fed i r n 



(21) 



At first sight, it would seem that each point r m would 
be within range of a large number of blip functions, so 
that many terms would have to be summed over for each 
r m in Eq. (^i"|). In fact, this is not so, provided the blip 
functions factorizc into Cartesian components in the way 
shown in Eq. (|l3|). To see this, assume that the blip grid 
and integration grid are cubic, let the blip-grid index £ 
correspond to the triplet (£ x ,£ y ,£ z ), and let the factor- 
ization of fe.(r) be written as: 



ft 



PiA x m) PtyiVm) PiA z m) , 



(22) 



where x m , y m and z m are the Cartesian components of 
r m (we suppress the index a for brevity) . The sum over I 
in Eq. (|2l| ) can then be performed as a sequence of three 
summations, the first of which is: 



( x ™) = h t x t y i z pe x (x m ) 



(23) 



The number of operations needed to calculate all 
these quantities 0£ y £ z (x m ) is just the number of points 
(£ x ,£ y , £ z ) on the blip grid times the number u- lrA of points 
x m for which pt x (x m ) is non-zero for a given £ x . This 
number u- lnt will generally be rather moderate, but the 
crucial point is that the number of operations involved is 
proportional only to v ln t and not to vf nt . Similar consid- 
erations will apply to the sums over £ y and £ z . 

It is worth remarking that since we have to calcu- 
late (pa (r m ) anyway, we have the option of calculating 
S a f3 by direct summation on the grid as well. In fact, 
T a (3 can also be treated this way, though here one must 
be more careful, since it is essential that its symmetry 
{T a (3 — Tf3a) be preserved by whatever scheme we use. 
This can be achieved by, for example, calculating the gra- 
dient V0 Q (r m ) analytically on the integration grid and 
then using integration by parts to express T a p as an in- 
tegral over V(p a ■ V0^. In the present work, we use the 
analytic forms for s a e,/3£> and t a e,pe'- 

A full linear-scaling calculation requires minimization 
of the total energy with respect to the quantities <p a and 
L a p- However, at present we are concerned solely with 
the representation of <p a , and the cut-off applied to L a p 
is irrelevant. For our practical tests of the blip- function 
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basis, we have therefore taken the L a p cut-off to infinity, 
which is equivalent to exact diagonalization of the Kohn- 
Sham equation. Apart from this, the procedure we use 
for determining the ground state, i.e. minimizing E with 
respect to the <f> a functions, is essentially the same as 
in our previous workEaEI We use conjugate-gradientsCJ 
minimization with respect to the blip-function coeffi- 
cients b a e- Expressions for the required derivatives are 
straightforward to derive using the methods outlined ear- 
lierM 



B. Practical tests 

We present here— ■ numerical tests both for the 
Appelbaum-HamanreJ local pseudopotential for Si used 
in Sec. 2 and for a standard Kerkera non-local pseudopo- 
tential for Si. The aims of the tests are: first, to show 
that the B-spline basis gives the accuracy in support- 
function calculations to be expected from our plane-wave 
calculations; and second to examine the convergence of 
E towards the exact plane-wave results as the region ra- 
dius i?reg is increased. For present purposes, it is not 
particularly relevant to perform the tests on large sys- 
tems. The tests have been done on perfect-crystal Si at 
the equilibrium lattice parameter, as in Sec. 2.3. 



h (A) 



E (eV/atom) 



0.4525 
0.3394 
0.2715 
0.2263 
0.1940 
0.1697 
0.1508 



-0.25565 
0.04880 
0.00818 
-0.01485 
0.00002 
-0.00270 
0.00000 



TABLE I. Total energy E as a function of integration grid 
spacing h using a region radius _R rcg = 2.715 A. The blip 
half-width &o was set to 0.905 A and the blip grid spacing 
used was 0.4525 A. The zero of energy is set equal to the 
result obtained with the finest grid. 



We have shown in Sec. 2.3 that a basis of B-splines 
having a half- width 6 — 1-0 A gives an error of ~ 
0.1 eV/atom if the blip spacing is ~ 0.45 A. For the 
present tests we have used the similar values &o = 0.905 A 
and a = 0.4525 A, in the expectation of getting this level 
of agreement with CASTEP plane-wave calculations in 
the limit R reg — > oo. To check the influence of integra- 
tion grid spacing h, we have made a set of calculations at 
different h using i? rcg = 2.715 A, which is large enough to 
be representative (see below). The results (Table 1) show 
that E converges rapidly with decreasing h, and ceases 
to vary for present purposes when h = 0.194 A. This 



confirms our expectation that h 



We have then 



used this grid spacing to study the variation of E with 
i? IO g, the results for which are given in Table 2, where we 
also compare with the plane-wave results. The extremely 
rapid convergence of E when i? reg exceeds ~ 3.2 A is very 
striking, and our results show that i? re g values yielding 
an accuracy of 10 -3 eV/atom are easily attainable. The 
close agreement with the plane- wave result fully confirms 
the effectiveness of the blip-function basis. As expected 
from the variational principle, E from the blip-function 
calculations in the R reg — > oo limit lies slightly above the 
plane-wave value, and the discrepancy of ~ 0.1 eV is of 
the size expected from the tests of Sec. 2.3 for (nearly) 
the present blip-function width and spacing. (We also re- 
mark in parenthesis that the absolute agreement between 
results obtained with two entirely different codes is useful 
evidence for the technical correctness of our codes.) 



i? reg (A) E (eV/atom) 

local pseudopotential non-local psuedopotential 



2.2625 
2.7150 
3.1675 
3.6200 
4.0725 



1.8659 
0.1554 
0.0559 
0.0558 
0.0558 



1.9653 
0.1507 
0.0396 
0.0396 
0.0396 



TABLE II. Convergence of the total energy E as a function 
of the region radius R IC g for silicon with a local and a non-local 
pseudopotential. The calculations were performed with a blip 
grid spacing of 0.4525 A and a blip half-width of 0.905 A in 
both cases. The zero of energy was taken to be the plane 
wave result obtained with each pseudopotential, with plane 
wave cutoffs of 250 and 200 eV respectively. 
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The results obtained in our very similar tests using 
the Kleinman-Bylander form of the Kerker pseudopoten- 
tial for Si are also shown in Table 2. In plane- wave cal- 
culations, the plane- wave cut-off needed for the Kerker 
potential to obtain a given accuracy is very similar to 
that needed in the Appclbaum-Hamann potential, and 
we have therefore used the same B-spline parameters. 
Tests on the integration-grid spacing show that we can 
use the value h = 0.226 A, which is close to what we have 
used with the local pseudopotential. The total energy 
converges in the same rapid manner for i? reg > 3.2 A, and 
the agreement of the converged result with the CASTEP 
value is also similar to what we saw with the Appclbaum- 
Hamann pseudopotential. 

IV. DISCUSSION 

In exploring the question of basis sets for linear-scaling 
calculations, we have laid great stress on the relation 
with plane-wave basis sets. One reason for doing this 
is that the plane-wave technique is the canonical method 
for pseudopotential calculations, and provides the easi- 
est way of generating definitive results by going to basis- 
set convergence. We have shown that within the linear- 
scaling approach the total energy can be taken to con- 
vergence by systematically reducing the width and spac- 
ing of a blip-function basis set, just as it can be taken 
to convergence by increasing the plane-wave cut-off in 
the canonical method. By analyzing the relation be- 
tween the plane- wave and blip- function bases, we have 
also given simple formulas for estimating the blip width 
and spacing needed to achieve the same accuracy as a 
given plane-wave cut-off. In addition, we have shown 
that the density of integration-grid points relates to the 
number of blip functions in the same way as it relates to 
the number of plane waves. Finally, we have seen that 
the blip-function basis provides a practical way of repre- 
senting support functions in linear-scaling calculations, 
and that the total energy converges to the plane-wave 
result as the region radius is increased. 

These results give useful insight into what can be ex- 
pected of linear-scaling DFT calculations. For large sys- 
tems, the plane-wave method requires a massive redun- 
dancy of information: it describes the space of occu- 
pied states using a number of variables of order N x M 
(N the number of occupied orbitals, M the number of 
plane waves) , whereas the number of variables in a linear- 
scaling method is only of order N x m (m the number of 
basis functions for each support function). This means 
that the linear-scaling method needs fewer variables than 
the plane- wave method by a factor m/M. But we have 
demonstrated that to achieve a given accuracy the num- 
ber of blip functions per unit volume is not much greater 
than the number of plane waves per unit volume. Then 
the factor m/M is roughly the ratio between the volume 
of a support region and the volume of the entire sys- 



tem. The support volume must clearly depend on the 
nature of the system. But for the Si system, we have 
seen that convergence is extremely rapid once the re- 
gion radius exceeds ~ 3.2 A, corresponding to a region 
volume of 137 A 3 , which is about 7 times greater than 
the volume per atom of 20 A 3 . In this example, then, 
the plane-wave method needs more variables than the 
linear-scaling method when the number of atoms N atom 
is greater than ~ 7, and for larger systems it needs more 
variables by a 'redundancy factor' of ~ N atom /7. (For 
a system of 700 atoms, e.g., the plane- wave redundancy 
factor would be ~ 100.) In this sense, plane-wave calcu- 
lations on large systems arc grossly inefficient. However, 
one should be aware that there are other factors in the 
situation, like the number of iterations needed to reach 
the ground state in the two methods. We are not yet in 
a position to say anything useful about this, but we plan 
to return to it. 

Finally, we note an interesting question. The impres- 
sive rate of convergence of ground-state energy with in- 
crease of region radius shown in Table 2 raises the ques- 
tion of what governs this convergence rate, and whether 
it will be found in other systems, including metals. We 
remark that this is not the same as the well-known ques- 
tion about the rate of decay of the density matrix p(r, r') 
as r — r'| — > oo, because in our formulation both <j) a (r) 
and L a fj play a role in the decay. Our intuition is that 
this decay is controlled by L a p. We hope soon to report 
results on the support functions for different materials, 
which will shed light on this. 
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APPENDIX: USING PLANE- WAVE 
CALCULATIONS TO TEST BLIP FUNCTIONS 

We explain here in more detail how a standard plane- 
wave code can be used to test the performance of blip- 
function basis sets (see Sec. 2.3). 

In a plane- wave calculation, the occupied orbitals if>i(r) 
are expressed as: 

Mr) = fi- 1/2 5> G exp(iG.r) , (Al) 

G 

where G goes over all reciprocal lattice vectors whose 
magnitude is less than the plane- wave cut-off G max , and 
ft is the volume per cell. (We assume T-point sampling 
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for simplicity.) The usual procedure is to minimize the 
total energy E with respect to the plane-wave coefficients 
c,Gi subject to the orthonormality constraints: 



E 



c iG c 3 G 



(A2) 



In a blip-function basis, the orbitals can be regarded 
as expanded in terms of the blip waves XG( r ) defined by: 



Xc(r) =Ag^2 ft(r) exp(iG • R £ ) , 



(A3) 



where the sum goes over all blip-grid points in one re- 
peating cell, and G is any of the reciprocal-lattice vec- 
tors of the repeating cell in the Brillouin zone associated 
with the blip grid. (The number of different G vectors 
labelling xg is thus the same as the number of points of 
the blip grid in one cell.) The xg(f) are orthogonal for 
different values of G, and we assume that the constant 
Ac, is chosen so that the XG( r ) are normalized: 



/ 



dr XG 1 XG 2 = <*Gi,G 2 • 



(A4) 



Then the orbitals are represented in terms of the blip 
waves as: 



(A5) 



G 



and the ground state is determined by minimizing E with 
respect to the blip coefficients biG subject to orthonor- 
mality: 



•jG 



Si-i ■ 



(A6) 



Now the qg can be expressed in terms of the biG by 
using the expansion of XG( r ) i n terms of plane waves: 

XG (r) =AGCo- 1 ^/(G + r)exp[ J (G + r)T] , (A7) 



different T. We can therefore reproduce the effect of a 
blip-function basis set by performing a standard plane- 
wave calculation in which the CjG+r are constrained to 
obey this relation. Note that in practice the plane-wave 
calculation must involve a plane- wave cut-off, so that the 
Ci G+r will actually be set to zero when |G + T| exceeds 
this cut-off. This is equivalent to setting /(G + T) equal 
to zero beyond the cut-off, and this amounts to a small 
smoothing of the blip function in real space. Allowing for 
this modification, the space spanned by the blip functions 
is a sub-space of the space spanned by the plane waves. 
The constraint on the Ci G+r is therefore a linear relation 
between these coefficients which expresses the fact that 
their component orthogonal to the sub-space vanishes, ft 
can be written as: 

- a G+r - f(G + T)J2 f* (G + T')d g+t' / 
r> 

]T|/(G + r')| 2 . (Aio) 

r' 

The modifications to the plane- wave code are therefore 
very simple. First, the initial guess for the orbitals must 
be chosen so as to satisfy this condition. Second, when 
the Ci G+r are updated on the iterative path to the self- 
consistent ground state, the changes of Cj G+r must be 
projected so as to satisfy this condition. This means that 
if the unconstrained calculation would have produced the 
changes <5ciG+r, these must be replaced by Scf G+r given 
by: 

6cj G+r = f(G + T)J2 f*{G + r')6c iG+r , I 
r' 

]T|/(G + r')| 2 . (Aii) 

r' 

This will produce exactly the same result as if we had 
worked directly with blip functions. 



where /(q) is the Fourier transform of /o(r), w is the vol- 
ume per blip-grid point, and T goes over the reciprocal- 
lattice vectors of the blip grid. From this, we have: 

^(r) = u- 1 A G b lG f(G + T) cxp [i(G + T) ■ r] , 
g r 

(A8) 



so that: 



Ci G+r = n^LJ^AabiGfiG + T) . 



(A9) 



Note that in interpreting this equation, G must be taken 
in the first Brillouin zone of the blip grid. 

This equation implies a relation between the plane- 
wave coefficients c^G+r for the same values of G but 
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